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In 2004 NASA began investigation of a robotic servicing mission for the Hubble Space 
Telescope (HST). Such a mission would require estimates of the HST attitude and rates 
in order to achieve a capture by the proposed Hubble robotic vehicle (HRV). HRV was to 
be equipped with vision-based sensors, capable of estimating the relative attitude between 
HST and HRV. The inertial HST attitude is derived from the measured relative attitude 
and the HRV computed inertial attitude. However, the relative rate between HST and 
HRV cannot be measured directly. Therefore, the HST rate with respect to inertial space 
is not known. Two approaches are developed to estimate the HST rates. Both methods 
utilize the measured relative attitude and the HRV inertial attitude and rates. First, a 
nonlinear estimator is developed. The nonlinear approach estimates the HST rate through 
an estimation of the inertial angular momentum. The development includes an analysis of 
the estimator stability given errors in the measured attitude. Second, a linearized approach 
is developed. The linearized approach is a pseudo-linear Kalman filter. Simulation test re- 
sults for both methods are given, including scenarios with erroneous measured attitudes. 

Even though the development began as an application for the HST robotic servicing mis- 
sion, the methods presented are applicable to any rendezvous/capture mission involving a 
non-cooperative target spacecraft. 

I. INTRODUCTION 

The Hubble Space Telescope (HST) was launched in 1990 and has undergone four servicing missions 

throughout its mission lifetime to replace instruments, sensors, solar arrays, power units, and cooling systems. 

‘Aerospace Engineer, julie.thienel@nasa.gov, phone:301-286-9033, and AIAA Member, 
t Professor of Aerospace Engineering, rmsanner@eng.umd.edu, phone:301-405-1928, and AIAA member. 
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Additional servicing will again be necessary to extend HST’s science, life. The batteries are predicted to fail 
as early as 2009, and the pointing control system may be reduced to a two-gyro mode by as early as 2007. 
On March 12, 2004 former NASA Administrator Sean O’Keefe asked the HST program to investigate robotic 
servicing of the HST to extend the science life. In addition to replacing the batteries and gyros, the servicing 
mission would install two new instruments and would provide the ability to safely de-orbit HST. 

The original robotic servicing mission concept included two vehicles, a de-orbit module (DM) and an 
ejection module (EM) with the robotic arm for servicing. The DM docks with HST and provides batteries 
and gyros for the stacked HST-DM configuration. The DM also provides the reentry capability for the 
stacked configuration. The EM carries the robotic arm for servicing, as well as the replacement instruments. 
Following the servicing, the EM separates from the HST-DM stack and safely de-orbits. An alternative 
concept to the robotic servicing mission was to provide a de-orbit only capability for HST. A modified DM 
docks with HST and provides for a safe re-entry of HST. In this work we refer to the vehicle docking with 
HST as the Hubble Robotic Vehicle (HRV). 

In order for the HRV to dock with HST, the HRV must match the rotation rates of HST. If either 
the batteries or gyros fail, HST will be tumbling with unknown rates. The HRV will be equipped with 
sensors capable of estimating the relative attitude between HST and HRV. The sensors utilize vision and 
feature recognition measurements to produce a relative quaternion. The HRV will also be equipped with 
star cameras and gyros which will provide the inertial orientation and body rates of the HRV. Combining 
the relative attitude and HRV inertial attitude provides a measurement of the HST attitude. However, the 
relative rates or HST body rates cannot be measured directly. 

We present and compare two approaches for estimating the HST body rates in the event that HST 
is tumbling and no a priori rate information is available. First a nonlinear estimator is developed. The 
nonlinear approach determines the HST rate through an estimation of the inertial angular momentum. The 
stability properties of the estimator are addressed first for perfect attitude measurements. Errors are then 
introduced into the HST measured attitude. The primary error source considered is the uncertainty in 
the relative attitude quaternion provided by the vision and feature recognition sensors. The results of the 
nonlinear approach are then compared to a pseudo-linear Kalman filter. 

Even though this work was originally developed for the HST/HRV scenario, it can easily be applied to 
other ‘non-cooperative’ rendezvous/capture scenarios. A ‘non-cooperative’ scenario is one in which the body 
rates of the target vehicle are not available, and the target vehicle does not carry any devices to improve the 
relative attitude estimation, such as reflectors or visible beacons. 

The paper is outlined as follows. The next section presents definitions and mathematical background. 
Section III presents the nonlinear estimation algorithm. Section IV discusses how measurement errors affect 
the nonlinear algorithm. Section V presents an overview of the pseudo-linear Kalman filter. Section VI 
includes results, followed by conclusions in the last section. 
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II. ATTITUDE AND ANGULAR RATE DEFINITIONS 


The attitude of a spacecraft can be represented by a quaternion, consisting of a rotation angle and unit 
rotation vector e, known as the Euler axis, and a rotation <j> about this axis so that 1 


Q = 


esin(^) 


£ 

CO 

O 

O 

1 


V 


(1) 


where q is the quaternion, partitioned into a vector part, e, and a scalar part, rj. The Hubble Space 
Telescope (HST) attitude quaternion is designated as q v , which defines the rotation from inertial to HST 
body coordinates. The Hubble Robotic Vehicle (HRV) attitude quaternion is designated as q h . 

The rotation, or attitude, matrix can be computed from the quaternion components as 1 

R - R(q) = (r ? 2 - £ T £)I 3 + ‘lee’ 1 ' - 2 t]S(e) (2) 

where I 3 is a 3x3 identity matrix and 5(e) is a matrix representation of the vector cross product operation. 
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Note also that R(q)e = e. 

A relative rotation between coordinate frames is computed as 2 
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Using the definition given in equation 3, the relative attitude quaternion from HRV body coordinates to 
HST body coordinates is then 


Qv ® Qh 


(4) 


The angular velocity of the HST body coordinates with respect to inertial space, resolved in HST body 
coordinates, is designated as u> v . Similarly the angular velocity of the HRV in HRV body coordinates is 
designated as u>h- The angular velocity of HST relative to HRV, given in HST body coordinates, is given as 


W” = - R r OJh 

where R r is the relative attitude matrix which transforms u>h into HST body coordinates. 


(5) 


III. HST NONLINEAR ANGULAR VELOCITY ESTIMATOR 

The following angular velocity estimator is intended for the scenario in which the HST batteries have died. 
No telemetry is available from HST. HRV is equipped with a quaternion star tracker, which gives q h , and 
some sensor system which produces the relative quaternion, q r . The approach does not yet account for any 
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errors in the measured quaternions, the HST inertia matrix, or the computation of the external torques on 
HST. The HST angular velocity is estimated in the inertial coordinate system through the estimation of the 
inertial angular momentum. The HST angular velocity in body coordinates is computed by a transformation 
of the inertial angular velocity. In the following developments, the time argument is omitted to simplify the 
notation. 

The angular momentum in inertial coordinates is defined as 


hi 


' T 1; ■ 


where 7i it) is the HST inertia matrix and u >i iV is the angular velocity, both in inertial coordinates. Note 
that Ij jV = R^I V R V , where R v is the HST attitude matrix defining the transformation from inertial to HST 
body coordinates. I v is the HST inertia matrix in body coordinates, assumed to be constant. The kinematic 
equation for the quaternion is given as 


~ 2 — 2^3(9u)Aj Rvhi,v 

where 


Q(q v ) = 

Vvh + S(£y) 


Qi(q v ) 
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-£ T 
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where, by inspection, Qi(q v ) = 77^,13 + S(e v ). The kinematic equation for the attitude matrix is 3 


( 6 ) 

(7) 
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Note that the q v (and therefore the attitude matrix R v ) is supplied by the HRV measured inertial attitude 
and the measured relative attitude. Using equation 4 


q v = q r ® qh 


(9) 


Similarly 

R/y R-yR-h 

Euler’s equation describes the dynamics of a rigid spacecraft. Euler’s equation for HST is given in inertial 
coordinates as 4 

hi, v = T iiV (10) 

Ti tV is the external torque acting on HST, resolved in inertial coordinates. 

Define the estimated angular momentum as 

hi t y — R^yUJi^y (H) 

where is the estimated angular velocity in inertial coordinates. The estimated angular velocity in body 
coordinates is then 

i*Jy RyUJi^y Ry I ^ y hi } y Iy Ryh^y 
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Using equation 1, define the predicted HST quaternion as 

Qv = 

Vv 

The attitude error is defined as the relative orientation between the predicted attitude q v and the attitude 
provided by the measured attitude, q v . The attitude error is 


Qv 


£v 

Vv 


= Qv®Q v 


( 12 ) 


Following the estimators proposed in Refs. 5, 6, and 4, a state estimator for the HST attitude and angular 
momentum is defined as 


h v = -Q{q v )R{q v (t)) T [I v 1 R v hi iV + ki v sign{fj v )] 


(13) 


h i>v = T i>v + —Ry I v e v sign(%) 


(14) 


The term R(q v (t)) T in equation 13 transforms the angular velocity terms from the body frame to the 
predicted attitude frame. The gain k is chosen as a positive constant. Similarly, the learning rate, a, is also 
a positive constant. Essentially, q v is a prediction of the attitude at time t, propagated with the kinematic 
equation using the estimated angular momentum. 

The kinematic equation for the attitude error quaternion, q v) has the same form as the quaternion 
kinematic equation in equation 6. The angular velocity associated with the attitude error quaternion is the 
difference between the angular velocity of the body coordinates and the angular velocity of the estimator 
coordinates (resolved in body coordinates). 7 Therefore, with equation 6, the definition given in equation 
7, equation 13, and noting that R(q v (t))i v = i v (since e v points along the eigenaxis of the rotation), the 
kinematic equation for q v is given as 


Qv 


1 


Qi(q v ) 


(/.y Ryhi^v ly Rv^i,V sign^^) 


Let hi yV = hi :V — hi jV . The derivative of h^ v is 




hi,v = -^Ryl v 1 e„sign(^) 


(15) 


(16) 


Note that the equilibrium states for 15 and 16 are 


hi. 


0 0 0 ±1 0 0 0 


Qv ,v i,v 

Proof of stability : Choose a Lyapunov function as (note that the time argument is omitted from the 
variables on the right side of the equation for clarity) 

1 - T . 1 \ {Vv - l) 2 + ele v ijv>0 

F(t) — 2^h'i l vh’i,v + 2 | 

\(fl v + l) 2 + ili v fj v <0 
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V(t ) is continuous. The derivative of V(t) is 


V(t) = -\ v h Uv + { 


(fjv - l)fjo + elio fj v >0 

(fj v + l)rj 0 + iyi 0 fjv < 0 


( 17 ) 


Noting that e^io + fj v f) 0 = 0 (including the left and right derivatives of fj v =0), for all t, equation 17 is 

F(t) = 

This establishes that hi tV , £ v , and fj v are globally, uniformly bounded. Moreover, V(t) is a continuous, twice 
differentiable function with 

■■ k 

V a (t) = --eZQ 1 (q v )[I~ 1 R v h i , v - ke v siga(rj v )} 

which is bounded. Barbalat’s lemma then shows that ||e„|| — > 0 as t —+ oo. s 

Since all signals in the estimator are bounded, the system 15 and 16 can be further analyzed, in the given 
format, as a linear time- varying system 9 x(t) = A(t)x(t ) where 


x(t) = 


£ v 

hi, v 


(18) 


A(t) = 


-%sign{fj v )Q 1 (q v ) \Qi{q v )I v l R v 
fsign 0 


where, by virtue of the above Lyapunov analysis, all terms in the matrix A{t) are known to be bounded 

, Theorem 4.5 and 


\fih 


0 


for all t > to- Rewriting V(t) as V(t) = —x(t) T C T Cx(t) < 0, where C = 
the discussion on pp. 626-628 in Ref. 9 shows that the equilibrium point x{t) = 0 of this equivalent system 
is exponentially stable if the pair (A(t),C) is uniformly completely observable (UCO). Since observability 
properties are unchanged under output feedback, 9 (^4(t), C) are UCO if the pair ( A(t ) —K ( t)C , C ) is uniformly 
observable for any piecewise, continuous and bounded matrix K(t). Choose K(t) as 


K(t) = 


fjv)Qi{q v ) 

V^sign(fi v )RT I~ l 


K(t) is piecewise continuous based on the following properties. Note from the above Lyapunov analysis that 
||e„|| — > 0. Since ||qj| 2 = 1 = [| || 2 + \fj v \ 2 for any time, t, there exists a time, T > 0, such that 1 1 77 ^, [ | > 0 
for all t > T. Since fj v therefore cannot pass through zero for t > T, sign(^) is constant for allt > T and, 
hence, K ( t ) is a piecewise continuous function of time. 

The state transition matrix for the pair {A(t) — K(t)C, C) is 


$(r,t) 


I 3 E(r,f) 

0 I 3 


( 19 ) 
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where E(r,f) = \ f t T Q 1 (q v (cr))I v 1 R v (a)da, with Qi(q v ) defined in equation 7. The observability Grammian 
is then 10 


pt+T pt+T 

W{t,t + T) = J $(r,t) T C T C${T,t)dT = Jt 


|l3 |S(r ,t) 

fS(r, t) T |E(r,t) r E(r,t) 


dr 


( 20 ) 


The system is UCO if there exists a T > 0 and positive constants a x < 00 , 0:2 > 0 such that, for all 
t > to ,o x I > W(t,t + T) > o 2 I. Using Lemma 13.4 of ref. 9, this property is assured if Qi(q v )I~ 1 R v and 
2 i(Qi(q v )Iy 1 R v ) are bounded, and there exist positive constants T 2 , Pi, and a finite constant /? 2 such that, 
for all t>to, 

pt+T2 

P 21 3>J t Rll- l Qi{q v ) T Qi{q v )I- x R v dT>pih (21) 

Qi(q v )I~ 1 R v is bounded, since all the signals in the estimator are bounded and the spacecraft inertia is 
bounded, hence the upper bound in 21 is satisfied. Substituting the equality Qi{q v ) T Qi{q v ) — I3 - 
into equation 21 results in 

pt+T 2 

00 > p 2 1 3 > / R t v I-\ I 3 - e v e T v )I~ l R v dT > p x l 3 > 0 (22) 


Recall that it has been shown that ||e„|| — > 0 asymptotically. Thus, for any 5 > 0, there exists a T X (S) > to 
such that ||£„|[ < 6 for all t > to + T x . Taking any <5 < 1, any T 2 >T X , and any 2 in M 3 

pt+T 2 

00 > (T 2 — Ti)A^ oa ,||z|| 2 > z T [ Jf (l-i v i T v )dr}z>(l-5 2 )(T 2 -T 1 )X 2 rnin \\zf>0 (23) 

where \ m ax and \ m in are the maximum and minimum eigenvalues of I~ x . Finally, jiQiiq^Iy 1 R v is 
bounded, since all the terms in 15 are bounded. This demonstrates the required UCO property. The 
PE condition is satisfied, and therefore i v and hi, v approach zero exponentially fast, i.e ||w„|| — > ||w„|| 
exponentially fast. □ 


IV. ERROR SOURCES 

The analysis provided in the previous section is for the ideal scenario with no errors in any of the 
measurements needed for the estimator. In reality, the measurements will be corrupted by errors. The 
primary error source considered here is the error in the relative attitude measurement used to derive the 
HST attitude quaternion. The resulting error in the HST attitude quaternion leads also to an error in the 
computation of the gravity gradient torque, the dominant external torque acting on HST. The following 
analysis considers how the HST attitude error affects the estimator stability and convergence properties. 
The measured HST quaternion, q v m , is written as 

Qv,m ^^err ® Qv (^4) 

where dq err is the error between the true HST quaternion, q v , and q ViTn ■ Recall that q v results from the 
product of the HRV relative attitude quaternion, q r , and the HRV inertial attitude quaternion, q h , shown 
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in equation 9. The error is written as 


SQerr 


From 24, the rotation matrices are given as 


£ er r 

7J err 


Rv,m — R(.Qv,m) R^Qerr R(& Q err) Rv 

Multiplying both sides of equation 24 by q v gives 


( 25 ) 


^Qerr ® Qv 


(26) 


Similarly, 

R(q v ,m) = R(Sq err )R(q v ) (27) 

Since the true quaternion is unknown, equations 13 and 14 cannot be implemented. Instead, q v m is used 
in place of q v , resulting in 


k v = 7;Q(Qv) R iQv,m( t )) T i I v lR v,mhi,v + fc£„, m sign('/7„ jln )] 


(28) 


A * * a O'* 1 ~ , / **• \ 

h-i } v = d" "2 Rv,m^v ^ v ,mSign(?7v ,m) 

Ti iV is the estimated external torque. Equations 25 and 27 are substituted into 28 

iu = \Q(qv)R{qv) T [R{ S Qerr) Tl v lR { 5 Q e rr) R vKv + kR{Sq err ) T £ Vtm S\gR{fl Vtm )\ 

Rearranging the terms results in 

Qv = ^QiQvjRiQvVtfv 1 Rvhi,v + ki v sign(fj v ) + A 1 h itV + A 2 ] 

where 

A! = [R(Sq err fl- l R(Sq err ) - I- l }R v 
A 2 = k[R(8q err ) T i Vim sign(fj Vi m) - e^sign^)] 


(29) 


(30) 


Ai and A 2 are both bounded. Ai results from the discrepancy in the HST inertia in the inertial frame due 
to the attitude measurement error. A 2 results directly from the attitude measurement error. 

Next, equation 29 is analyzed. Recall that the true angular momentum is driven by the external torques 
acting on HST, shown in equation 10. The dominant external torque acting on HST is gravity gradient 
torque. The gravity gradient torque in inertial coordinates is given as 11 


T itV = ^-S{rl v )R T v I v R v rl v 


(31) 


where p is the earth’s gravitational parameter, Ti tV is the magnitude of the HST inertial position vector, and 
is the inertial position unit vector. The estimated torque in equation 29 is computed as 


Ti,v = ^S(rl v )Rl m I v R Vtm rl v 


( 32 ) 
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Substituting equation 25 into 32, fi yV is written as 


where 


Ti, v = 4£ S{rl v )R T v I v R v rl v + A 3 = T<,„ + A 3 


3/t 


A 3 = ^fS(rl v )[Rl m I v R v , m - Rll v R v ]rl v 

^i,v 

Substituting 33 and 25 into 29 results in 


where 


h i>v = T i>v + A 3 + -R* v I v 1 e„sign(^ w ) + A 4 


A 4 = -^Rl[R{Sq e rr) T I v ^u.mSign (fjv.m) ~ Iy ^Sigll^,,)] 


Combining equations 30 and 34 with 6 and 10, we get 


Qv,m 


\ Q(q v (t))(I~ 1 R v h i:V - k£ v sign(fj v )) + A[h itV + A^ 



-f .R^J-^sign^) + A' 3 + A; 


where A[ = -%Q(q v (t)) A x , A 2 = -§Q(q v (t)) A 2 , A' 3 = -A 3 , and A' 4 = -A 4 
Equation 35 has the form 

x(t) = f(t,x(t)) +d(x,t) 

where f(x,t) has been shown to be exponentially stable and 


d(t) = 


= di(t) + d 2 (t) 


di (t) = 


( 33 ) 


(34) 


(35) 


Ai hi iV + A 2 
A'+A^ 

where 

r t . 

A^i^ 

0 

and d 2 (t) contains the remaining terms, which are all bounded. Since d 2 (t) is bounded, the exponentially 
stable system is robust to this perturbation. The concern, however, is with d\{t) which contains part of the 
state and hence cannot be assumed bounded a priori. If exponential stability is preserved with c?i(f), then 
the system will be robust to the combined disturbance d(t). 

The nominal system is exponentially stable. Therefore, according to the Converse Lyapunov Theorem, 
a Lyapunov function and positive constants ci, C 2 , c 3 , and c 4 exist for the nominal system and satisfy the 
following 9 

ci||x(f)|| < V p {t) < c 2 ||x(f)|| 

V P (t) < -c 3 \\x(t)f 

< C 4 ||x(t)|| 


dV p (t) 


dx 
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Using V p (t ) as the Lyapunov function candidate for the system perturbed by di(t), the derivative of V p (t) 
satisfies 

F P (t) < - C3 ||*(t)ll 2 + 11^11 ll<*i(*)ll 

which becomes 

V P (t) < -c 3 ||*(f)|| 2 + c 4 ||as(t)||J|rfi(i)|| (36) 

Note that we can write 


<f 


^1 ^i,v) 


^1 ^i,v d” ^i,v 

0 


0 


0 


A[hi tV can be moved into d 2 (f ) since the angular momentum of HST will be bounded. Equation 36 is then 


T P (i)<- C3 |kWI| 2 + C4||*(t)|||| 


-A [h itV 
0 


V p (t) < c 3 1| cc(t) || 2 + c 4 J||x(t)|| 2 (37) 

where || || <6. If <5 is sufficiently small, exponential stability is preserved, and the system is robust to the 

combined perturbation d(t) = di(t) + d 2 (t). Thus, the proposed algorithm is still stable in the presence of 
attitude measurement errors, provided that those errors are sufficiently small. 

Adding d 2 (t) back into the Lyapunov analysis, equation 37 becomes 

V p (t) < -c 3 ||a:(t)|| 2 + c 4 <5||x(t)|| 2 + c 4 |jcc(i) || || rf 2 (*) II (38) 


The bound on || || , which establishes the size of 6 , and the bound on ||cZ 2 (i) || are defined next. From 

equation 2, R(6q err ) can be written as 


^(^q e rr) — ( Verr ^err^err) I 3 T 2£ err £ err c 2j] err S{ K S err ) — I 3 2£ err £ err 2 £ err £ err 27] err S{s err ) 


= Is + Y (, 1 


erri & err 


(39) 


and 

||^(^/err5 ^err)|| = 2||£errj| 

Also, note that ||Q(g„(f))|| = 1, ||-R„|| = ||-R(gJ|| = 1, ||e„|| < 1, ||e„ im || < 1, \\I V \\ < A max , and || /“^l < 
The || || is then bounded as 

II || < ^HAil < HJ|££!zJL(i + ||e err .||) < S 

The ||d 2 (t)|| is bounded as 


ll*(*)ll < II a 2 II + 1 ! A 3 1 | + HA4II + ||Ai||||fc i ,„|| mas 


(40) 
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with 


IIAall < fc||e e rr|| 

HA^H < 12MA T 3 al — (1 + lkTrll) 

' i,v 

II A4H < 7— ~ ll £ err|| 

^ min 

The angular momentum of HST is bounded as 


— ll^6,vll — ^max [1^6 ,^, max || 


Extensive simulations of HST predict the maximum angular velocity as approximately 0.004 rad/sec. 12 
Therefore, 

11 a' mu. 11 ^ 0.008A mO 3,||£ er . r || / . i M _ 11 ■, 

II Ai j| || hi )V || rnax 51 , (IT ll^err j| J 

^ min 

Applying Young’s inequality to equation 38, and combining terms gives 


V P (t) < -(ic 3 - C4<5)il*(t)|| 2 + ^WMt)\\ 2 
Substituting all the bounds into equation 40, equation 41 is written as 

V p (t) < ~(^c 3 - c A 5)\\x(t)\\ 2 + -^\\£ err \\ 2 (a + b\\e e 


(41) 


(42) 


where a -k + j- 2 - + ■- 2 ^J - a - £ + P-°° 8A ™ a:c and b - ^ 


.008A m acc 


. Recall that ||e e rr|| = sin|0 err . 


Assuming the <p err is small, ||£ e 


\4>err- Equation 42 is then 


Vp{t) < -(i C3 - C4<5)||*(f)|| 2 + §-/err(a + \b<j> err ) 2 


Neglecting powers of <t> e rr greater than 2 gives 


,1 


2„2 


crc; 


Vp{t) < -(|c 3 - c 4< 5)||x(t)|| 2 + 

Computing the time average of equation 44, and taking the limit as T —* oo, is 


(43) 


(44) 


lim sup 

T — >oo 


i f T 

tL w 


a 2 c\ 


t)|| dr < limsup . 

r-.<x> 8 c 3 (c 3 - C 4 S) J 0 


/ <perr(r) 2 dT 

Jo 


If the error process is ergodic, the ensemble average is equivalent to the time average. The RMS bound is 
given as 

\\hi,v\\RMS < II^COII-RMS < c<7 (45) 


where u is the standard deviation of the error process, and all the constants are combined together as 
r _ .act 

8C3(C3-C45) ' 

Equation 45 serves to demonstrate that the angular momentum estimator is bounded, and converges 
to a ball sized by the bound given in 45 (given that the attitude is corrupted by errors meeting the above 
statistical assumptions) . Equation 45 sizes the errors only to within the original unknown constants of the 
converse Lyapunov analysis. The characteristics of the estimator error will be addressed through numerical 
simulations. 
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In the case that the effects of attitude measurement errors cannot be assumed small (that is 5 in eqn 
(37) is not sufficiently small), the potential instability resulting from di(t) can be removed through a minor 
modification of the estimator. This is accomplished by introducing a ‘leakage’ term in equation 29. 


where 


hi^y — 7 z , v ~b ^ rn ^ v ^{hi )V ^fli )V 


&{h'i y v') — \ 


(J 0 if ||hj )t ,j[ || Umax ||h.j jV || 


(46) 


0 otherwise 

where the scalar a Q > 0 and \\hi^ v \\ max is the upper bound on the HST angular momentum. This modification 
serves to force the magnitude of the angular momentum estimate to decrease if ever it exceeds a known prior 
bound on the feasible magnitude of the actual HST momentum. 

The additional term in 46 introduces an additional term into the analysis of 38 


-t ? 




Applying Young’s inequality results in 

__ rp ^ 1 ~ 1 

crh i<v h i<v < --ff\\h itV \\ 2 + -o-llhi^U 2 

Since hi V is assumed to be bounded, the second term acts as another bounded disturbance, which can be 
combined with d, 2 . The first term, however, acts to stabilize the possible high noise instability identified 
above. In particular, if <r 0 > 2a5 then the system again exponentially converges to a ball about the origin, 
regardless of the size of <5. The size of the ball to which the state converges can be bounded similarly to the 
development above. 

Although it looks like this technique adds an additional disturbance, possibly stabilizing the system but 
adding extra error sources, in fact this is not the case. Expanding the added term instead as 


— rr(A^ jV hi jV ) hi iV — rrhi^ v hi f y ~ cr\\ h i,v|| 


Or 


Recall from above the a > 0 only if ||/ii,„|| > ||fci,u||- Thus the added term is always negative; essentially the 
extra constraint in the estimation law for the angular momentum can only help convergence. Note that if 
the size of the perturbation 5 due to the measurement errors is sufficiently large to require use of leakage to 
stabilize the algorithm, the resulting RMS prediction errors are likely to be quite poor. Since the leakage can 
only improve the convergence, however, it is a useful safety feature to include in the algorithm, regardless of 
the noise level. 
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V. HST PSEUDO-LINEAR ANGULAR VELOCITY ESTIMATOR 


The pseudo-linear Kalman filter is based on the model presented in reference 13. A brief overview of the 
filter models is given first, followed by a summary of the actual algorithm equations. 

The state of the filter is 

X: r qv 

where q v is the HST inertial attitude and uj v is the HST angular velocity in HST body coordinates. The 
quaternion is propagated according to equation 6 as 

1 


Qv = j Q(Qv)“v 

The angular velocity is propagated using Euler’s equation as 

^ v = I v [S(I v u: v )u> v + T& )V ] 

where T b]V is the gravity gradient torque from equation 32, transformed to body coordinates as 

T b<v = ^S(r% v )I v R v rZ v 

r i,v 

where r£ v = R v r^ v . Equation 49 can be written in a pseudo- linear form as 

T b , v = ^S(rl v )I v M(rZ v ,q v )q v = F gg q v 


(47) 


(48) 


(49) 


(50) 


where M(rf v , q v ) is defined in reference 14. Equation 50 is substituted into 48, and equations 47 and 48 are 
then augmented as 

X = F{X)X + w(t) (51) 

where w(t), a zero-mean white noise process, is added to account for the lack of knowledge in the terms 
given in equations 47 through 50. F(X) is defined as 


F(X) = 


0 


\Q(.Qv) 


Iv l F gg I-'Sih «„) 


(52) 


The measurement used in the filter is q v m given in equation 24. The measurement model for a given 
time tk is simply 


L 0 


Xk = HXk + Vk 


Qv,m,k 

where I 4 is a 4x4 identity matrix and Vk is a zero-mean white measurement noise. 
The filter algorithm is implemented discretely as follows: 

Propagation: 


X k+1 {-) = $ k X k (+) 
= e F(X k (+))At 
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At — tjc+1 

Pk+i(-) = *kPk{+)*I + Q 

where P is the filter covariance, Q is the process noise covariance matrix, and (— ) and (+) indicate a priori 
and a posteriori values, respectively. 

Update: 


X k (+) = A fc (-) + K k (q Vtmtk - q k (-)) 

K k = HP k {-){HP k {-)H T + R)- 1 
P k (+) = (I - K k H)P k (-)( I - K k H) T + K k RKl 

where K k is the filter gain matrix at time, t k , and R is the measurement noise covariance matrix. 

VI. SIMULATION RESULTS 

The algorithms outlined in the previous sections are tested with a simulation. The simulation includes 
two orbits of data, based on an actual HST ephemeris. The HST inertia is 12 


36046 

-706 

1491 

-706 

86868 

449 

1491 

449 

93848 


Both algorithms are initially tested without any errors. In both cases, the initial attitude quaternions 
are identity, q v = q v = [0,0,0, 1). The initial HST angular velocity estimate is zero, Cj v = [0,0,0], and the 
true initial angular velocity is = [-0.04,-0.01,0.14] deg/sec. The gains are chosen as k = 0.005 and 
a = 9e 5 . The noise matrices in the pseudo-linear filter are R = la;10 _e l4 and Q = lxlO -13 !?, and the initial 
covariance is P t o = I7. 

Figures 1(a) and 1(b) are examples of the resulting angular velocity error from the nonlinear estimator 
and the pseudo-linear filter, respectively. The pseudo-linear algorithm converges very quickly, but the error 
oscillates. The nonlinear estimator has a longer convergence time, however the errors after convergence are 
very small. Next, the measured attitude is chosen randomly with a 15 degree uncertainty. Figures 2(a) 
and 2(b) are samples of the angular velocity error for the nonlinear estimator and the pseudo-linear filter, 
respectively, using the same erroneous measured attitudes. As with the clean example above, the pseudo- 
linear filter converges faster, but is noisier than the nonlinear estimator in steady state. Table 1 lists the 
average of the final RMS angular velocity error for each of the 100 test cases, for both the nonlinear estimator 
and the pseudo-linear Kalman filter. The average errors are smaller for the nonlinear estimator on all three 
axes. 

Next, the leakage term of equation 46 is added to the nonlinear estimator, with cr 0 = 1. In order to see 
the improvement in the estimation, the initial angular velocity estimate is increased such that the initial 
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Figure 1. Angular Velocity Errors, No Attitude Measurement Error, (a) Nonlinear, (b) Pseudo-Linear 


angular momentum estimate exceeds the maximum anticipated angular momentum. Figure 3(a) shows the 
angular velocity errors with the ‘leakage’ term, and 3(b) shows the same scenario without the leakage term. 
In both cases the nonlinear estimator converges, however the additional term does improve the convergence 
time. Figure 3(c) shows the same scenario with the pseudo-linear filter. The pseudo-linear filter does not 
perform as well with very large initial errors. 

Finally, figure 4 shows the RMS angular momentum at the end of two orbits for each of the 100 test 
cases. These results could be utilized in analyzing the numerical bound expressed in equation 45. 
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Figure 2. Angular Velocity Errors, Random 15 Degree Measurement Error, (a) Nonlinear, (b) Pseudo-Linear 


VII. CONCLUSIONS 

A nonlinear algorithm is developed to estimate the rotation rates for a non-cooperative target vehicle. The 
nonlinear algorithm determines the rotation rate through an estimation of the inertial angular momentum. 
The algorithm, designed for the Hubble Robotic Servicing mission, applies in particular to the scenario in 
which the batteries have died and HST is tumbling. The HRV design includes vision and feature recognition 
sensors capable of producing a relative attitude quaternion. Combining the relative attitude with the HRV 
inertial attitude produces the measured HST attitude for the estimation algorithm. The nonlinear algorithm 
is compared to a pseudo-linear Kalman filter. 

The two rate estimation approaches are compared using two orbits of HST data. Initially, the algorithms 
are compared with perfect measured attitudes. Both algorithms estimate the HST rotation rate. The 
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Table 1. Average Final RMS Angular Velocity Errors 



Nonlinear (deg/sec) 

Pseudo-Linear (deg/sec) 

&x,RMS 

0.00164 

0.00516 

tdy,RMS 

0.00164 

0.00558 

&z,RMS 

0.00127 

0.00555 


pseudo-linear algorithm converges quickly, but oscillates slightly in steady state. The nonlinear approach 
is slower to converge; the steady state errors, however, do not oscillate and are very nearly zero. The two 
approaches are then compared for 100 different cases, with different random measured attitudes in each case. 
Again, both approaches estimate the rotation rate quite well. However, the nonlinear approach produces 
consistently smaller final rotation rate errors, especially given large initial errors in the estimates. 

Future work will focus on improving the accuracy of the relative attitude quaternion computed from 
the vision and feature recognition sensors. Additional error sources, such as uncertainties in the spacecraft 
inertia, should be addressed, as well as a more rigorous study of the nonlinear estimator gains. Also, the 
algorithms will be coupled with a control scheme to study the closed loop behavior during a rendezvous or 
docking scenario. 
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Figure 3. Angular Velocity Errors, Large Initial Angular Velocity Estimate, (a) Nonlinear with Leakage, (b) 
Nonlinear without Leakage (c) Pseudo-Linear 
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Figure 4. Final Angular Momentum Error RMS 
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